##############################################
# Waves of securitisation: the rise, fall and resurgence of citizenship stripping regulations in Europe, 1960 - 2022
# Luuk van der Baaren, Maria Gerdes and Maarten Vink
# Statelessness and Citizenship Review
# Code to reproduce plots in article and Annex
# Code prepared by Maarten Vink

#start with clean workspace
rm(list=ls(all=TRUE))

#load packages
library(scales)
library(dplyr)
library(tidyr)
library(ggplot2)
library(gmodels)
library(grid)
library(gridtext)
library(ggrepel)
library(janitor)
library(readxl)
library(viridis)
library(openxlsx)
library(data.table)
library(writexl)
library(forcats)
library(ggpubr)


# read in data
d03 <- read_excel("Waves_of_securitization_data_final.xlsx", sheet = "L03")
d04 <- read_excel("Waves_of_securitization_data_final.xlsx", sheet = "L04")
d07 <- read_excel("Waves_of_securitization_data_final.xlsx", sheet = "L07")
d08 <- read_excel("Waves_of_securitization_data_final.xlsx", sheet = "L08")
d07_allegiance <- read_excel("Waves_of_securitization_data_final.xlsx", sheet = "L07_allegiance")
d07_security <- read_excel("Waves_of_securitization_data_final.xlsx", sheet = "L07_security")
d07_terrorism <- read_excel("Waves_of_securitization_data_final.xlsx", sheet = "L07_terrorism")
d07_keywords <- read_excel("Waves_of_securitization_data_final.xlsx", sheet = "L07_keywords")


# Descriptive plots L03 - L08 regulation

# Import data and selected only coded columns and 
# Remove countries not independent and code unknown legislation as NA
# Bind all data frames and order data into a single file
# Note some columns only for L07 so for other modes these are treated as NA's
dta <- bind_rows(d03, d04, d07, d08) |>
  filter(value_bin != 95) |>
  filter(value_cat != 95) |>
  mutate(value_bin = ifelse(value_bin == 98, NA, value_bin)) |>
  mutate(value_cat = ifelse(value_cat == 98, NA, value_cat)) |>
  mutate(intro_cat = as.factor(intro_cat)) |>
  arrange(country, year, mode_id)
View(dta)
summary(dta) #8879 obs of 18 vars
# iso3               iso2              region            country               year        mode_id           mode_desc        
# Length:8879        Length:8879        Length:8879        Length:8879        Min.   :1960   Length:8879        Length:8879       
# Class :character   Class :character   Class :character   Class :character   1st Qu.:1978   Class :character   Class :character  
# Mode  :character   Mode  :character   Mode  :character   Mode  :character   Median :1996   Mode  :character   Mode  :character  
# Mean   :1994                                        
# 3rd Qu.:2009                                        
# Max.   :2022                                        
# 
# article           procedure           value_bin        value_cat      specification       condition        
# Length:8879        Length:8879        Min.   :0.0000   Min.   :0.0000   Length:8879        Length:8879       
# Class :character   Class :character   1st Qu.:0.0000   1st Qu.:0.0000   Class :character   Class :character  
# Mode  :character   Mode  :character   Median :0.0000   Median :0.0000   Mode  :character   Mode  :character  
# Mean   :0.3091   Mean   :0.4302                                        
# 3rd Qu.:1.0000   3rd Qu.:1.0000                                        
# Max.   :1.0000   Max.   :2.0000                                        
# NA's   :32       NA's   :32                                            
#                                         intro_cat      intro_bin       value_bin1      value_bin2    stateless_prot 
# Provision introduced                   :  19           Min.   :0.000   Min.   :0.000   Min.   :0.000   Min.   :0.000  
# Scope extension                        :  11           1st Qu.:0.000   1st Qu.:0.000   1st Qu.:0.000   1st Qu.:0.000  
# Provision abolished                    :   6           Median :0.000   Median :0.000   Median :0.000   Median :0.000  
# Provision introduced - state succession:   6           Mean   :0.024   Mean   :0.023   Mean   :0.048   Mean   :0.422  
# Scope limitation                       :   4           3rd Qu.:0.000   3rd Qu.:0.000   3rd Qu.:0.000   3rd Qu.:1.000  
# (Other)                                :   3           Max.   :1.000   Max.   :1.000   Max.   :2.000   Max.   :1.000  
# NA's                                   :8830           NA's   :6659    NA's   :6659    NA's   :6659    NA's   :6659  
head(dta)

#remove NA's and selected only coded columns
dt <- dta |>
  filter(!is.na(value_bin)) |>
  select(c(2,4,5,6,10,11,13,14,15,18))
summary(dt) #8847 obs of 10 vars without NA's
# iso2             country               year        mode_id            value_bin        value_cat       condition        
# Length:8847        Length:8847        Min.   :1960   Length:8847        Min.   :0.0000   Min.   :0.0000   Length:8847       
# Class :character   Class :character   1st Qu.:1978   Class :character   1st Qu.:0.0000   1st Qu.:0.0000   Class :character  
# Mode  :character   Mode  :character   Median :1996   Mode  :character   Median :0.0000   Median :0.0000   Mode  :character  
#                                       Mean   :1994                      Mean   :0.3091   Mean   :0.4302                     
#                                       3rd Qu.:2009                      3rd Qu.:1.0000   3rd Qu.:1.0000                     
#                                       Max.   :2022                      Max.   :1.0000   Max.   :2.0000                     
# 
#                                     intro_cat      intro_bin     stateless_prot 
# Provision introduced                   :  19   Min.   :0.000   Min.   :0.000  
# Scope extension                        :  11   1st Qu.:0.000   1st Qu.:0.000  
# Provision abolished                    :   6   Median :0.000   Median :0.000  
# Provision introduced - state succession:   6   Mean   :0.024   Mean   :0.421  
# Scope limitation                       :   4   3rd Qu.:0.000   3rd Qu.:1.000  
# (Other)                                :   3   Max.   :1.000   Max.   :1.00
# NA's                                   :8798   NA's   :6632    NA's   :6632 
head(dt)

# rename columns
data <- dt |>
  select(c(1,2,3,4,5,6,8,10)) |>
  tidyr::pivot_wider(
    names_from  = c(mode_id), # Can accommodate more variables, if needed.
    values_from = c(value_bin, value_cat, intro_cat, stateless_prot)) |>
  rename("L03_bin" = "value_bin_L03")|>
  rename("L04_bin" = "value_bin_L04")|>
  rename("L07_bin" = "value_bin_L07")|>
  rename("L07_intro" = "intro_cat_L07") |>
  rename("L07_stat_prot" = "stateless_prot_L07") |>
  rename("L08_bin" = "value_bin_L08")|>
  rename("L03_cat" = "value_cat_L03")|>
  rename("L04_cat" = "value_cat_L04")|>
  rename("L07_cat" = "value_cat_L07")|>
  rename("L08_cat" = "value_cat_L08") |>
  select(-c(12,13,15,16,17,19)) # select only relevant columns
View(data)


# Figure 1: data visualisation showing number or % of countries by mode
# with one or more of these loss/deprivation grounds
# select 3 years to visualise trend

#descriptives
data |>
  filter(year==1960 | year==1990 | year==2022) |>
  count(year)
# year     n
# 1  1960    28
# 2  1990    30
# 3  2022    42

#combined overall measure
#combined measures
#combined L03 or L04 or L07 or L08
data$L03_4_7_8<- ifelse(data$L03_bin == 1 | data$L04_bin == 1 | data$L07_bin == 1 | data$L08_bin == 1, "Yes",
                    "No")
data |>
  filter(year==1960 | year==1990 | year==2022) |>
  count(year, L03_4_7_8)  |>
  group_by(year) |>
  mutate(freq = percent(n / sum(n))) 
# year L03_4_7_8     n freq 
# 1  1960 No           12 43%  
# 2  1960 Yes          15 54%  
# 3  1960 NA            1 4%    # Due to missing for Austria, L03/L04 in 1960-64, FR: 1960
# 4  1990 No           10 33%  
# 5  1990 Yes          20 67%  
# 6  2022 No           12 29%  
# 7  2022 Yes          30 71%  

#change into long format & remove _bin
#count by mode
loss_count_bin3 <- data |> 
  select(c("iso2", "year", "L03_bin", "L04_bin", "L07_bin", "L08_bin")) |> 
  pivot_longer(names_to = "mode", values_to = "category", L03_bin:L08_bin)
loss_count_bin3$mode = sub("\\_bin.*","", as.character(loss_count_bin3$mode))
View(loss_count_bin3) #8860 obs of 4 variables

#recode values 0 and 1 in whole dataframe
class(loss_count_bin3$category)
loss_count_bin3$category <- as.character(loss_count_bin3$category)
loss_count_bin3$category[loss_count_bin3$category=='0'] <- 'No'
loss_count_bin3$category[loss_count_bin3$category=='1'] <- 'Yes'
loss_count_bin3$mode[loss_count_bin3$mode=='L03'] <- 'Service in foreign army'
loss_count_bin3$mode[loss_count_bin3$mode=='L04'] <- 'Other service to a foreign country'
loss_count_bin3$mode[loss_count_bin3$mode=='L07'] <- 'Disloyalty or state security'
loss_count_bin3$mode[loss_count_bin3$mode=='L08'] <- 'Other offences'
View(loss_count_bin3) #8860 obs of 4 variables


# Data for Fig 1
Fig1_loss_count_bin_mode_yr <- loss_count_bin3 |> 
  na.omit() |>
  count(year, mode, category) |>
  group_by(mode, year) |> 
  mutate(percent = percent(n/sum(n),1))
Fig1_loss_count_bin_mode_yr


# Figure 1
plot_modes_loss_count_bin_nt <- Fig1_loss_count_bin_mode_yr |>
  na.omit() |>
  filter(year==1960 | year==1990 | year==2022) |> 
  ggplot(aes(x=factor(mode), y = n,fill=factor(category)))+geom_bar(stat="identity")+
  coord_flip()+
  facet_wrap(~year)+
  geom_text(aes(label = percent), size = 5, position=position_stack(vjust=0.5))+ 
  theme(legend.position='none')+
  theme_minimal(base_size = 22)+
  theme(strip.text.x = element_text(size = 20, face = "bold", hjust = 0, vjust = -1))+
  scale_fill_manual(name = "", 
                    values = c("#8fccc6", "#c1cdcd"))+
  theme(legend.position='top', legend.justification='left')+
  guides(fill=guide_legend("Citizenship deprivation ground present in country", reverse=TRUE))+
  theme(axis.text.x = element_text(size = 20))+
  theme(axis.text.y = element_text(size = 20))+
  xlab("")+
  ylab("\nnumber of countries (and %)\n")+
  labs(caption = '')
plot_modes_loss_count_bin_nt

jpeg('Fig1.modes.loss.count.selyrs.jpeg',  width = 18, height = 12, units = 'in', res = 200)
plot_modes_loss_count_bin_nt
dev.off()


  # Figure 2: visualisation of trends by mode

# calculate the mean share of countries with a deprivation rule, by year, by mode
daty <- dt |> 
  # na.omit() |>
  group_by(year, mode_id) |>
  summarise(mean_yr = mean(value_bin, na.rm = TRUE)) |>
  ungroup()
daty

# adjust labels
daty$mode_id[daty$mode_id=='L03'] <- 'Service in foreign army'
daty$mode_id[daty$mode_id=='L04'] <- 'Other service to a foreign country'
daty$mode_id[daty$mode_id=='L07'] <- 'Disloyalty or state security'
daty$mode_id[daty$mode_id=='L08'] <- 'Other offences'


############ Figure 2. L03-L08 trend 

jpeg('Fig2.loss.trend.jpeg',  width = 12, height = 8, units = 'in', res = 200)
ggplot(daty, aes(fill=mode_id, y=mean_yr, x=year)) + 
  geom_bar(position="dodge", stat="identity")+
  facet_wrap(~mode_id)+
  theme_minimal(base_size = 22)+
  theme(legend.position="none") +
  xlab("")+
  ylab("")+
  scale_fill_manual(name = "", 
                    values = c("#8fccc6", "#76a7a2","#5e837f","#47615e"))+
  scale_x_continuous(breaks = seq(1960, 2022,10))+
  scale_y_continuous(labels = scales::percent, breaks = seq(0,1.0,0.2))+
  labs(title="", 
       subtitle='Share of countries with citizenship deprivation ground',
       caption = '')
dev.off()

# Robustness check: Figure A1
# facet by old/new country
# check trend with only countries that exist from 1960 to 2023
dt1 <- dt |>
  group_by(country) |>
  arrange(year) |>
  mutate(newcountry = ifelse(min(year)<1990, "Countries independent before 1990", "Countries independent after 1989"))
# recode Czech Republic as successor of Czechoslovakia
dt1$newcountry[dt1$country =="Czech Republic"] <- "Countries independent before 1990"
# recode Russia as successor of Soviet Union
dt1$newcountry[dt1$country =="Russia"] <- "Countries independent before 1990"

dt1 <- dt1 |>
  mutate(newcountry = factor(newcountry, levels = c("Countries independent before 1990", "Countries independent after 1989")))

dat1y <- dt1 |> 
  # na.omit() |>
  group_by(year, newcountry, mode_id) |>
  summarise(mean_yr = mean(value_bin, na.rm = TRUE)) |>
  ungroup()
dat1y

# labels
dat1y$mode_id[dat1y$mode_id=='L03'] <- 'Service in foreign army'
dat1y$mode_id[dat1y$mode_id=='L04'] <- 'Other service to a foreign country'
dat1y$mode_id[dat1y$mode_id=='L07'] <- 'Disloyalty or treason'
dat1y$mode_id[dat1y$mode_id=='L08'] <- 'Other offences'

jpeg('FigA1.loss.trend.oldnew.jpeg',  width = 10, height = 13, units = 'in', res = 200)
ggplot(dat1y, aes(fill=mode_id, y=mean_yr, x=year)) + 
  geom_bar(position="dodge", stat="identity")+
  facet_grid(mode_id~newcountry)+
  theme_minimal(base_size = 16)+
  theme(legend.position="none") +
  xlab("")+
  ylab("")+
  scale_fill_manual(name = "", 
                    values = c("#8fccc6", "#76a7a2","#5e837f","#47615e"))+
  scale_x_continuous(breaks = seq(1960, 2022,10))+
  scale_y_continuous(labels = scales::percent, breaks = seq(0,1.0,0.2))+
  labs(title="", 
       subtitle='',
       caption = '')
dev.off()


# Fig 3: count number of security-related loss modes per country

# count loss modes per country, select columns with loss modes
loss_count_bin <- data |> 
  select(c("iso2","country",  "year", "L03_bin", "L04_bin", "L07_bin", "L08_bin")) |> 
  group_by(iso2, country, year) |>
  mutate(n.loss.modes = rowSums (across(everything())))
View(loss_count_bin) #2215 obs of 8 variables
summary(loss_count_bin)


# Figure 3: Plot number of loss modes by country
plot_data <- loss_count_bin |>
  na.omit() |>
  filter(year==1960 | year==1990 | year==2022) |> 
  group_by(year, n.loss.modes) |> 
  tally() |> 
  mutate(percent = percent(n/sum(n),1))
View(plot_data)


# Plot 
jpeg('Fig3.loss.modes.count.jpeg',  width = 27, height = 8, units = 'in', res = 200)
plot_data |>
  ggplot(aes(x = n.loss.modes, y = n, fill = factor(n.loss.modes))) +
  geom_bar(stat = "identity") +
  facet_wrap(~year)+
  geom_text(aes(label = percent), size = 5, vjust = 0, hjust = -0.05)+ 
  coord_flip()+
  theme(legend.position='none')+
  theme_minimal(base_size = 20)+
  theme(strip.text.x = element_text(size = 20, face = "bold", hjust = 0, vjust = -1))+
  theme(legend.position='none')+
  theme(axis.text.x = element_text(size = 20))+
  theme(axis.text.y = element_text(size = 20))+
  scale_fill_manual(name = "", 
                    values = c("#8fccc6", "#76a7a2","#5e837f","#47615e", "#30413f"))+
  xlab("")+
  ylab("\nnumber of countries (and %)")+
  labs(subtitle = 'Number of citizenship deprivation grounds by country') 
dev.off()


# Figure 4: changes dotplot L07

# new dataframe with changes only
data_changes <- data |> filter(L07_intro != "No substantive change") |>
  mutate(L07_intro_rec = L07_intro) |> # add new variable so we can recode it
  mutate(year = as.numeric(year)-1) |> # deduct 1 year in order to plot in the correct year of change (and not from the subsequent year, as in overall GLOBALCIT approach)
  select(iso2, country, year, L07_intro, L07_intro_rec, L07_cat)

# recode state succession changes in order to include these also
data_changes$L07_intro_rec[data_changes$L07_intro_rec =="Provision abolished - state succession"] <- "Provision abolished"
data_changes$L07_intro_rec[data_changes$L07_intro_rec =="Provision introduced - state succession"] <- "Provision introduced"
#order type of changes
data_changes$L07_intro_rec <- factor(data_changes$L07_intro_rec,levels = c("Provision abolished", "Scope limitation", 
                                                      "Scope extension", "Provision introduced"))

# Plot and save
jpeg('Fig4_L07.changes.type.jpeg',  width = 10, height = 5, units = 'in', res = 200)
data_changes |>
  mutate(change_type = ifelse(L07_intro_rec == 'Provision abolished' | L07_intro_rec == 'Scope limitation', 'Statelessness risk decreases',
                              ifelse(L07_intro_rec == 'Provision introduced' | L07_intro_rec == 'Scope extension', 'Statelessness risk increases', NA))) |>
  ggplot(aes(year, fill = L07_intro_rec))+
  geom_dotplot(binwidth=1, stackgroups = TRUE, binpositions = "all")+
  facet_wrap(~change_type, ncol = 1)+
  scale_y_continuous("",
                     limits=c(0,4), breaks = seq(0,4,2))+
  scale_x_continuous("",
                     limits=c(1960, 2022), breaks = seq(1960, 2022, 10))+
  scale_fill_manual("Type of changes:", values = c("#00703c", "#88d8c0", "#bda9a9",  "#fa7b62")) +
  theme_bw(base_size = 16)+
  theme(legend.position="top")+
  coord_fixed(ratio=1)+
  labs(title="", 
       caption = '')
dev.off()



# Fig 5: categorical: discriminatory provisions
#change into long format & remove _bin
#count by mode
loss_count_cat <- data |> 
  select(c("iso2", "year", "L03_cat", "L04_cat", "L07_cat", "L08_cat")) |> 
  pivot_longer(names_to = "mode", values_to = "category", L03_cat:L08_cat)
loss_count_cat$mode = sub("\\_cat.*","", as.character(loss_count_cat$mode))
View(loss_count_cat) #8860 obs of 4 variables

#recode values 0 and 1 in whole dataframe
class(loss_count_cat$category)
loss_count_cat$category <- as.character(loss_count_cat$category)
loss_count_cat$category[loss_count_cat$category=='0'] <- 'No provision'
loss_count_cat$category[loss_count_cat$category=='1'] <- 'General provision'
loss_count_cat$category[loss_count_cat$category=='2'] <- 'Only certain categories of citizens'
loss_count_cat$mode[loss_count_cat$mode=='L03'] <- 'Service in foreign army'
loss_count_cat$mode[loss_count_cat$mode=='L04'] <- 'Other service to a foreign country'
loss_count_cat$mode[loss_count_cat$mode=='L07'] <- 'Disloyalty or state security'
loss_count_cat$mode[loss_count_cat$mode=='L08'] <- 'Other offences'
View(loss_count_cat) #8860 obs on 4 variables


# Data for Fig 5
Fig5_loss_count_cat_mode_yr <- loss_count_cat |> 
  na.omit() |>
  count(year, mode, category) |>
  group_by(mode, year) |>
  mutate(percent = percent(n/sum(n),1))
Fig5_loss_count_cat_mode_yr

Fig5_loss_count_cat_mode_yr$category <- as.factor(Fig5_loss_count_cat_mode_yr$category)
class(Fig5_loss_count_cat_mode_yr$category)


# Plot and save
jpeg('Fig5.modes.loss.cat.count.jpeg',  width = 20, height = 8, units = 'in', res = 200)
Fig5_loss_count_cat_mode_yr |>
  na.omit() |>
  filter(year==1960 | year==2022) |> 
  mutate(category = fct_relevel(category, 
                                "No provision", "General provision", "Only certain categories of citizens")) |>
  ggplot(aes(x=factor(mode), y = n,fill=factor(category)))+geom_bar(stat="identity")+
  coord_flip()+
  facet_wrap(~year)+
  geom_text(aes(label = percent), size = 6, position=position_stack(vjust=0.5))+ 
  theme(legend.position='none')+
  theme_minimal(base_size = 22)+
  theme(strip.text.x = element_text(size = 22, face = "bold", hjust = 0, vjust = -1))+
  scale_fill_manual(name = "", 
                    values = c("#c1cdcd", "#83bfb9","#a8e6df"))+
  theme(legend.position='top', legend.justification='left')+
  guides(fill=guide_legend("Deprivation ground", reverse=TRUE))+
  theme(axis.text.x = element_text(size = 22))+
  theme(axis.text.y = element_text(size = 22))+
  xlab("")+
  ylab("\nnumber of countries (and %)\n")+
  labs(caption = '')
dev.off()

# Fig 6: Statelessness risk? 

# make new statelessness protection variable (categorical)
data <- data |>
  mutate(L07_stat_prot2 = ifelse(L07_bin == 0, "No provision",
                                 ifelse(L07_bin == 1 & L07_stat_prot == 0, 'Provision with statelessness risk',
                                        ifelse(L07_bin == 1 & L07_stat_prot == 1, 'Provision with protection against statelessness',NA))))
data$L07_stat_prot2 <- factor(data$L07_stat_prot2,levels = c("No provision", "Provision with statelessness risk", 
                                                             "Provision with protection against statelessness"))
levels(data$L07_stat_prot2)

Fig6_L07_count_stat_prot <- data |> 
  group_by(year) |>
  count(L07_stat_prot2) |>
  mutate(percent = percent(n/sum(n),1))
Fig6_L07_count_stat_prot

# Plot and save
jpeg('Fig6.L07.trend.statelessness.prot.jpeg',  width = 14, height = 10, units = 'in', res = 200)
data |> 
  group_by(year) |>
  count(L07_stat_prot2) |>
  mutate(percent = n/sum(n)) |>
  ggplot(aes(year, percent, fill = L07_stat_prot2)) +
  geom_area(alpha=0.6 , linewidth=.5, colour="white")+
  scale_y_continuous(labels = scales::percent)+
  scale_x_continuous("",
                     limits=c(1960,2025), breaks = seq(1960, 2025, 20))+
  scale_fill_manual(name = "", 
                    labels = c("No provision", "Provision with statelessness risk", "Provision with protection against statelessness"),
                    values = c("#6ba7a1", "#30413f", "#807532"))+
  ylab("") +
  xlab("")+
  theme_minimal()+
  theme(text = element_text(size = 22))+
  theme(legend.position="top")+
  labs(title="Citizenship deprivation based on disloyalty or treason", 
       subtitle='Share of countries with type of provision\n',
       caption = '') +
  theme(plot.title = element_text(hjust = 0, size = 24)) +
  theme(plot.subtitle = element_text(hjust = 0, size = 18))
dev.off()


### Figure: prevalence of terms

# Figure 7a: Prevalence of disloyalty-related terms 

## prepare data 
dat_all <- d07_allegiance |>
  group_by(year) |>
  filter(procedure != "n.a.") |>
  count(keywords_allegiance_bin) |>
  mutate(percent = n/sum(n))

# df with only 1's
dat_all1 <- dat_all |>
  filter(keywords_allegiance_bin == 1)
dat_all1

# Figure 7: Prevalence of disloyalty-related terms 
Fig7a_allegiance <- dat_all1 |> 
  ggplot(aes(year, percent, fill=factor(keywords_allegiance_bin))) +
  geom_area(alpha=0.6 , linewidth=.5, colour="white")+
  scale_y_continuous(labels = scales::percent,
                     limits=c(0,0.95))+
  scale_x_continuous("",
                     limits=c(1960,2025), breaks = seq(1960, 2025, 20))+
  ylab("") +
  xlab("")+
  theme_minimal()+
  scale_fill_manual(name = "", values = c("#6ba7a1"))+
  theme(text = element_text(size = 22))+
  theme(legend.position="none")+
  labs(title="", 
       subtitle='Allegiance\n',
       caption = '')
Fig7a_allegiance

# Figure 7a: Prevalence of disloyalty-related terms 

## prepare data 
dat_sec <- d07_security |>
  group_by(year) |>
  filter(procedure != "n.a.") |>
  count(keywords_security_bin) |>
  mutate(percent = n/sum(n)) 

# dataframe with only 1's
dat_sec1 <- dat_sec|>
  filter(keywords_security_bin == 1)
dat_sec1

# Figure 7b: Prevalence of security-related terms 
Fig7b_security <- dat_sec1 |> 
  ggplot(aes(year, percent, fill=factor(keywords_security_bin))) +
  geom_area(alpha=0.6 , linewidth=.5, colour="white")+
  scale_y_continuous(labels = scales::percent,
                     limits=c(0,0.95))+
  scale_x_continuous("",
                     limits=c(1960,2025), breaks = seq(1960, 2025, 20))+
  scale_fill_manual(name = "", values = c("#7c9d99"))+
  ylab("") +
  xlab("")+
  theme_minimal()+
  theme(text = element_text(size = 22))+
  theme(legend.position="none")+
  labs(title="", 
       subtitle='Security\n',
       caption = '')
Fig7b_security

# Figure 7c: Prevalence of terrorism-related terms 

## prepare data 
dat_ter <- d07_terrorism |>
  group_by(year) |>
  filter(procedure != "n.a.") |>
  count(keywords_terrorism_bin) |>
  mutate(percent = n/sum(n))

dat_ter1 <- dat_ter |>
  filter(keywords_terrorism_bin == 1)
dat_ter1

# Figure 7: Prevalence of disloyalty-related terms 
Fig7c_terrorism <- dat_ter1 |> 
  ggplot(aes(year, percent, fill=factor(keywords_terrorism_bin))) +
  geom_area(alpha=0.6 , linewidth=.5, colour="white")+
  scale_y_continuous(labels = scales::percent,
                     limits=c(0,0.95))+
  scale_x_continuous("",
                     limits=c(1960,2025), breaks = seq(1960, 2025, 20))+
  scale_fill_manual(name = "", values = c("#899391"))+
  ylab("") +
  xlab("")+
  theme_minimal()+
  theme(text = element_text(size = 22))+
  theme(legend.position="none")+
  labs(title="", 
       subtitle='Terrorism\n',
       caption = '')
Fig7c_terrorism

# combi plot
library(patchwork)
jpeg('Fig7_keywords.jpeg',  width = 14, height = 10, units = 'in', res = 200)
(Fig7a_allegiance + Fig7b_security +Fig7c_terrorism) +
  plot_annotation(title = "The language of deprivation provisions regarding disloyalty and state security",
                  subtitle = "Share of citizenship laws with relevant provisions containing related terms") & 
  theme(plot.title = element_text(hjust = 0, size = 24)) &
  theme(plot.subtitle = element_text(hjust = 0, size = 18))
dev.off()



# Save data used for plots
# write multiple sheets
dataset_names <- list('Fig1_loss_count_bin_mode_yr' = Fig1_loss_count_bin_mode_yr,
                      'Fig2_loss_trend' = daty,
                      'Fig3_loss_count' = plot_data,
                      'Fig4_changes' = data_changes,
                      'Fig5_discr' = Fig5_loss_count_cat_mode_yr,
                      'Fig6_statl' = Fig6_L07_count_stat_prot,
                      'Fig7a_allegiance' = dat_all,
                      'Fig7b_security' = dat_sec,
                      'Fig7c_terrorism' = dat_ter)
write.xlsx(dataset_names, file = 'Waves_of_secur_plots_data.xlsx')

### END
